use mergedmerged, clear

encode state, gen(id)
sort id year
xtset id year

/*interpolating missing data*/
ipolate blackpct year, by(state) gen(blacki)
ipolate pvi year, by(state) gen(pvii)
ipolate urban year, by(state) gen(urbani)

/*standardizing variables*/
egen blackis=std(blacki)
egen pviis=std(pvii)
egen urbanis=std(urbani)
egen bipcs=std(bipc)
egen dems=std(demo)
egen reps=std(repu)
egen citis=std(citi6)
egen libs=std(liberal)
egen cons=std(conserv)
egen urb40s=std(urb40)
replace south=L.south if missing(south)
replace north=L.north if missing(north)
replace ncent=L.ncent if missing(ncent)
replace west=L.west if missing(west)
gen bipt1=bipt+1
gen urb40sl=L30.urb40s
gen citisq=citis^2
gen pvisq=pviis^2

/*summary statistics*/
sort year
by year: sum brifle blocks bballist basslt bmags bgunban citis
/*nonvarying policies excluded from gsem*/

/*table 2*/
gsem (Guns -> bminage bmult bssp botheft bdesign bdtr bconst btrain, probit) (Guns -> bmg bss bsbs bsbr baow bgunban bwait bdealer bpriv blic bregis bretent bipcs boci bcci, regress) (Guns -> bipt1, gamma) if year==1986, latent(Guns) startvalues(iv)
est store A
gsem (Guns -> bminage bmult bssp botheft bdesign bconst btrain, probit) (Guns -> brifle bballist blocks basslt bmags bdtr bmg bss bsbs bsbr baow bwait bdealer bpriv blic bregis bretent bipcs boci bcci, regress) (Guns -> bipt1, gamma) if year==2016, latent(Guns) startvalues(iv)
est store B
est table A B, se stats(bic)

/*figure 1*/
local i=1986
forvalues i=1986/1988 {
	quietly gsem (Guns -> bminage bmult bssp botheft bdesign bdtr bconst btrain, probit) (Guns -> bmg bss bsbs bsbr baow bgunban bwait bdealer bpriv blic bregis bretent bipcs boci bcci, regress) (Guns -> bipt1, gamma) if year==`i', latent(Guns) startvalues(iv)
    predict g`i' if e(sample), latent(Guns) se(gs`i')
}

local i=1989
forvalues i=1989/1993 {
	quietly gsem (Guns -> bminage bmult bssp botheft bdesign bconst btrain, probit) (Guns -> basslt bmags bdtr bmg bss bsbs bsbr baow bgunban bwait bdealer bpriv blic bregis bretent bipcs boci bcci, regress) (Guns -> bipt1, gamma) if year==`i', latent(Guns) startvalues(iv)
    predict g`i' if e(sample), latent(Guns) se(gs`i')
}

local i=1994
forvalues i=1994/1999 {
	quietly gsem (Guns -> bminage bmult bssp botheft bdesign bconst btrain, probit) (Guns -> blocks basslt bmags bdtr bmg bss bsbs bsbr baow bgunban bwait bdealer bpriv blic bregis bretent bipcs boci bcci, regress) (Guns -> bipt1, gamma) if year==`i', latent(Guns) startvalues(iv)
    predict g`i' if e(sample), latent(Guns) se(gs`i')
}

local i=2000
forvalues i=2000/2002 {
	quietly gsem (Guns -> bminage bmult bssp botheft bdesign bconst btrain, probit) (Guns -> bballist blocks basslt bmags bdtr bmg bss bsbs bsbr baow bgunban bwait bdealer bpriv blic bregis bretent bipcs boci bcci, regress) (Guns -> bipt1, gamma) if year==`i', latent(Guns) startvalues(iv)
    predict g`i' if e(sample), latent(Guns) se(gs`i')
}

local i=2003
forvalues i=2003/2013 {
	quietly gsem (Guns -> bminage bmult bssp botheft bdesign bconst btrain, probit) (Guns -> brifle bballist blocks basslt bmags bdtr bmg bss bsbs bsbr baow bgunban bwait bdealer bpriv blic bregis bretent bipcs boci bcci, regress) (Guns -> bipt1, gamma) if year==`i', latent(Guns) startvalues(iv)
    predict g`i' if e(sample), latent(Guns) se(gs`i')
}

local i=2014
forvalues i=2014/2016 {
	quietly gsem (Guns -> bminage bmult bssp botheft bdesign bconst btrain, probit) (Guns -> brifle bballist blocks basslt bmags bdtr bmg bss bsbs bsbr baow bwait bdealer bpriv blic bregis bretent bipcs boci bcci, regress) (Guns -> bipt1, gamma) if year==`i', latent(Guns) startvalues(iv)
    predict g`i' if e(sample), latent(Guns) se(gs`i')
}

gen gindex=.
gen ghi=.
gen glo=.

local i=1986
forvalues i=1986/2016 {
	replace gindex=g`i' if year==`i'
	replace ghi=g`i'+gs`i'*1.96 if year==`i'
	replace glo=g`i'-gs`i'*1.96 if year==`i'
}

graph box gindex if year>1985 & year<2017, over(year, label(angle(90) labsize(*0.8))) ytitle("Gun Control Index")

/*figure 2*/
twoway rcap ghi glo pviis if year==1986 || scatter gindex pviis if year==1986, mlabel(abb) mlabsize(vsmall) mlabpos(0) msymbol(none) xtitle("Democratic + Green Lean (PVI)") ytitle("Gun Control Index") || lfit gindex pviis if year==1986
/*figure 3*/
twoway rcap ghi glo pviis if year==2016 || scatter gindex pviis if year==2016, mlabel(abb) mlabsize(vsmall) mlabpos(0) msymbol(none) xtitle("Democratic + Green Lean (PVI)") ytitle("Gun Control Index") || lfit gindex pviis if year==2016

save gindex, replace

/*figure 4*/
local i=1986
forvalues i=1986/2014 {
	corr gindex median if year==`i'
	scalar rho`i'=r(rho)
}
matrix define rhos=(`=rho1986' \ `=rho1987' \ `=rho1988' \ `=rho1989' \ `=rho1990' \ `=rho1991' \ `=rho1992' \ `=rho1993' \ `=rho1994' \ `=rho1995' \ `=rho1996' \ `=rho1997' \ `=rho1998' \ `=rho1999' \ `=rho2000' \ `=rho2001' \ `=rho2002' \ `=rho2003' \ `=rho2004' \ `=rho2005' \ `=rho2006' \ `=rho2007' \ `=rho2008' \ `=rho2009' \ `=rho2010' \ `=rho2011' \ `=rho2012' \ `=rho2013' \ `=rho2014')
preserve
drop if year<1986
drop if year>2014
drop if state~="Alabama"
mkmat year, mat(Y)
clear
matrix define cor=(rhos,Y)
matrix colnames cor = rho year 
svmat double cor, names(col)
twoway bar rho year , ytitle("Correlation Coefficient")
restore

/*test correlations identical*/

preserve
drop if year~=1986
save first, replace
restore
drop if year~=2014
save second, replace
append using first
gen d=0
replace d=1 if year==2014
reg gindex median c.median#d, robust

/*/*policy lib loadings*/
gsem (Liberalism -> abortion_m~d death_pen drugs_mari~o education_~a environmen~l environme~as gambling_l~n genderrigh~n genderrigh~s immigrati~ge income_taxes labor_prev~s labor_righ~k labor_stat~e v86 regulatio~es regulati~lts r~otorcycl~s w_ec_access w_gayright~n bmult bconst , probit) (Liberalism -> z_labor_un~o x_sales_ta~s w_educatio~e w_abor~1983_ bdtr bgunban bwait bpriv blic bregis , regress) (Guns -> bminage bmult bssp botheft bdesign bconst btrain, probit) (Guns -> bdtr bmg bss bsbs bsbr baow bgunban bwait bdealer bpriv blic bregis bretent bipcs boci bcci, regress) (Guns -> bipt1, gamma), latent(Guns Liberalism) startvalues(iv) group(d) ginvariant(cons coef errvar scale means loading)
*/

clear
use gindex
sort id year
xtset id year


/*opinion liberalism & partisanship*/
gsem (Guns -> bminage bmult bssp botheft bdesign bdtr bconst btrain, probit) (Guns -> bmg bss bsbs bsbr baow bgunban bwait bdealer bpriv blic bregis bretent bipcs boci bcci, regress) (Guns -> bipt1, gamma) (blackis urbanis citis citisq south -> Guns, ) if year==1986, latent(Guns) startvalues(iv)
matrix b=e(b)
est store A

gsem (Guns -> bminage bmult bssp botheft bdesign bdtr bconst btrain, probit) (Guns -> bmg bss bsbs bsbr baow bgunban bwait bdealer bpriv blic bregis bretent bipcs boci bcci, regress) (Guns -> bipt1, gamma) (blackis urbanis pviis pvisq south -> Guns, ) if year==1986, latent(Guns) from(b, skip)
est store B

gsem (CCW -> btrain, probit) (CCW -> bipcs boci bcci, regress) (CCW -> bipt1, gamma) (blackis urbanis citis south -> CCW, ) if year==1986, latent(CCW) startvalues(iv)
est store Z
predict CCW86 if e(sample), latent(CCW)
matrix b=e(b)

gsem (CCW -> btrain, probit) (CCW -> bipcs boci bcci, regress) (CCW -> bipt1, gamma) (blackis urbanis pviis south -> CCW, ) if year==1986, latent(CCW) from(b, skip)
est store Y

gsem (Guns -> bminage bballist bmult bssp botheft bdesign bdtr bconst btrain, probit) (Guns -> bmg bss bsbs bsbr baow bwait bdealer bpriv blic bregis bretent brifle blocks bballist basslt bmags bipcs boci bcci, regress) (Guns -> bipt1, gamma) (L6.blackis L6.urbanis L3.citis south -> Guns) if year==2016, latent(Guns) startvalues(iv) 
matrix b=e(b)
est store C

gsem (Guns -> bminage bballist bmult bssp botheft bdesign bdtr bconst btrain, probit) (Guns -> bmg bss bsbs bsbr baow bwait bdealer bpriv blic bregis bretent brifle blocks bballist basslt bmags bipcs boci bcci, regress) (Guns -> bipt1, gamma) (L6.blackis L6.urbanis pviis south -> Guns) if year==2016, latent(Guns) from(b, skip)
est store D

gsem (CCW -> btrain, probit) (CCW -> bipcs boci bcci, regress) (CCW -> bipt1, gamma) (L6.blackis L6.urbanis L3.citis south -> CCW, ) if year==2016, latent(CCW) startvalues(iv) 
matrix b=e(b)
est store X
predict CCW16 if e(sample), latent(CCW)

gsem (CCW -> btrain, probit) (CCW -> bipcs boci bcci, regress) (CCW -> bipt1, gamma) (L6.blackis L6.urbanis pviis south -> CCW, ) if year==2016, latent(CCW) from(b, skip)
est store W

est table A B C D, se stats(bic)
/*est table Z Y X W, se stats(bic)*/

/*table 3*/
gen CCW=CCW86
replace CCW=CCW16 if missing(CCW)
drop CCW86 CCW16

reg CCW L30.CCW L30.gindex L30.citis if year==2016
hettest
reg CCW L30.CCW L30.gindex L30.citis if year==2016, robust
est store M
reg CCW L30.CCW L30.gindex L30.pviis if year==2016, robust
est store N
reg CCW L30.CCW L30.gindex L30.citis L3.citis L30.pviis pviis if year==2016, robust
est store O
est table M N O, se stats(bic)

drop g19* g20* gs19* gs20* gindex ghi glo

/*figure 5*/
replace citis=L.citis if missing(citis)
replace citisq=L.citisq if missing(citisq)
replace urbanis=L.urbanis if missing(urbanis)
replace blackis=L.blackis if missing(blackis)
local i=1986
forvalues i=1986/1988 {
	quietly gsem (Guns -> bminage bmult bssp botheft bdesign bdtr bconst btrain, probit) (Guns -> bmg bss bsbs bsbr baow bgunban bwait bdealer bpriv blic bregis bretent bipcs boci bcci, regress) (Guns -> bipt1, gamma) (blackis urbanis citis citisq south -> Guns, ) if year==`i', latent(Guns) startvalues(iv)
    predict g`i' if e(sample), latent(Guns) se(gs`i')
    est store E`i'
	lincom citis+citisq
	matrix b`i'=e(b)
	matrix r`i'=b`i'[1,49]
	matrix q`i'=b`i'[1,50]
	matrix c`i'=b`i'[1,51]
	scalar c`i'=c`i'[1,1]
	scalar r`i'=r`i'[1,1]
	scalar q`i'=q`i'[1,1]
	matrix d`i'=b`i'[1,52]
	scalar d`i'=d`i'[1,1]
	matrix s`i'=e(V)
	matrix t`i'=s`i'[51,51]
	scalar t`i'=t`i'[1,1]
	scalar c`i'hi=`=c`i''+1.96*sqrt(`=t`i'')
	scalar c`i'lo=`=c`i''-1.96*sqrt(`=t`i'')
	matrix u`i'=s`i'[52,52]
	scalar u`i'=u`i'[1,1]
	scalar d`i'hi=`=d`i''+1.96*sqrt(`=u`i'')
	scalar d`i'lo=`=d`i''-1.96*sqrt(`=u`i'')
	matrix v`i'=s`i'[49,49]
	scalar v`i'=v`i'[1,1]
	matrix w`i'=s`i'[50,50]
	scalar w`i'=w`i'[1,1]
	scalar r`i'hi=`=r`i''+1.96*sqrt(`=v`i'')
	scalar r`i'lo=`=r`i''-1.96*sqrt(`=v`i'')
	scalar q`i'hi=`=q`i''+1.96*sqrt(`=w`i'')
	scalar q`i'lo=`=q`i''-1.96*sqrt(`=w`i'')
	nlcom _b[Guns:citis]-2*_b[Guns:citisq], post
	matrix liblo`i'=e(b)
	matrix lls`i'=e(V)
	scalar liblo`i'lo=liblo`i'[1,1]-1.96*sqrt(`=lls`i'[1,1]')
	scalar liblo`i'hi=liblo`i'[1,1]+1.96*sqrt(`=lls`i'[1,1]')
	est restore E`i'
	nlcom _b[Guns:citis]+2*_b[Guns:citisq], post
	matrix libhi`i'=e(b)
	matrix lhs`i'=e(V)
	scalar libhi`i'lo=libhi`i'[1,1]-1.96*sqrt(`=lhs`i'[1,1]')
	scalar libhi`i'hi=libhi`i'[1,1]+1.96*sqrt(`=lhs`i'[1,1]')
}

local i=1989
forvalues i=1989/1993 {
	quietly gsem (Guns -> bminage bmult bssp botheft bdesign bconst btrain, probit) (Guns -> basslt bmags bdtr bmg bss bsbs bsbr baow bgunban bwait bdealer bpriv blic bregis bretent bipcs boci bcci, regress) (Guns -> bipt1, gamma) (blackis urbanis citis citisq south -> Guns, ) if year==`i', latent(Guns) startvalues(iv)
    predict g`i' if e(sample), latent(Guns) se(gs`i')
    est store E`i'
	lincom citis+citisq
    matrix b`i'=e(b)
	matrix c`i'=b`i'[1,55]
	scalar c`i'=c`i'[1,1]
	matrix d`i'=b`i'[1,56]
	scalar d`i'=d`i'[1,1]
	matrix r`i'=b`i'[1,53]
	matrix q`i'=b`i'[1,54]
	scalar r`i'=r`i'[1,1]
	scalar q`i'=q`i'[1,1]
	matrix s`i'=e(V)
	matrix t`i'=s`i'[55,55]
	scalar t`i'=t`i'[1,1]
	scalar c`i'hi=`=c`i''+1.96*sqrt(`=t`i'')
	scalar c`i'lo=`=c`i''-1.96*sqrt(`=t`i'')
	matrix u`i'=s`i'[56,56]
	scalar u`i'=u`i'[1,1]
	scalar d`i'hi=`=d`i''+1.96*sqrt(`=u`i'')
	scalar d`i'lo=`=d`i''-1.96*sqrt(`=u`i'')
	matrix v`i'=s`i'[53,53]
	scalar v`i'=v`i'[1,1]
	matrix w`i'=s`i'[54,54]
	scalar w`i'=w`i'[1,1]
	scalar r`i'hi=`=r`i''+1.96*sqrt(`=v`i'')
	scalar r`i'lo=`=r`i''-1.96*sqrt(`=v`i'')
	scalar q`i'hi=`=q`i''+1.96*sqrt(`=w`i'')
	scalar q`i'lo=`=q`i''-1.96*sqrt(`=w`i'')
	nlcom _b[Guns:citis]-2*_b[Guns:citisq], post
	matrix liblo`i'=e(b)
	matrix lls`i'=e(V)
	scalar liblo`i'lo=liblo`i'[1,1]-1.96*sqrt(`=lls`i'[1,1]')
	scalar liblo`i'hi=liblo`i'[1,1]+1.96*sqrt(`=lls`i'[1,1]')
	est restore E`i'
	nlcom _b[Guns:citis]+2*_b[Guns:citisq], post
	matrix libhi`i'=e(b)
	matrix lhs`i'=e(V)
	scalar libhi`i'lo=libhi`i'[1,1]-1.96*sqrt(`=lhs`i'[1,1]')
	scalar libhi`i'hi=libhi`i'[1,1]+1.96*sqrt(`=lhs`i'[1,1]')
}

local i=1994
forvalues i=1994/1999 {
	quietly gsem (Guns -> bminage bmult bssp botheft bdesign bconst btrain, probit) (Guns -> blocks basslt bmags bdtr bmg bss bsbs bsbr baow bgunban bwait bdealer bpriv blic bregis bretent bipcs boci bcci, regress) (Guns -> bipt1, gamma) (blackis urbanis citis citisq south -> Guns, ) if year==`i', latent(Guns) startvalues(iv)
    predict g`i' if e(sample), latent(Guns) se(gs`i')
    est store E`i'
	lincom citis+citisq
    matrix b`i'=e(b)
	matrix c`i'=b`i'[1,57]
	scalar c`i'=c`i'[1,1]
	matrix d`i'=b`i'[1,58]
	scalar d`i'=d`i'[1,1]
	matrix r`i'=b`i'[1,55]
	matrix q`i'=b`i'[1,56]
	scalar r`i'=r`i'[1,1]
	scalar q`i'=q`i'[1,1]
	matrix s`i'=e(V)
	matrix t`i'=s`i'[57,57]
	scalar t`i'=t`i'[1,1]
	scalar c`i'hi=`=c`i''+1.96*sqrt(`=t`i'')
	scalar c`i'lo=`=c`i''-1.96*sqrt(`=t`i'')
	matrix u`i'=s`i'[58,58]
	scalar u`i'=u`i'[1,1]
	scalar d`i'hi=`=d`i''+1.96*sqrt(`=u`i'')
	scalar d`i'lo=`=d`i''-1.96*sqrt(`=u`i'')
	matrix v`i'=s`i'[55,55]
	scalar v`i'=v`i'[1,1]
	matrix w`i'=s`i'[56,56]
	scalar w`i'=w`i'[1,1]
	scalar r`i'hi=`=r`i''+1.96*sqrt(`=v`i'')
	scalar r`i'lo=`=r`i''-1.96*sqrt(`=v`i'')
	scalar q`i'hi=`=q`i''+1.96*sqrt(`=w`i'')
	scalar q`i'lo=`=q`i''-1.96*sqrt(`=w`i'')
	nlcom _b[Guns:citis]-2*_b[Guns:citisq], post
	matrix liblo`i'=e(b)
	matrix lls`i'=e(V)
	scalar liblo`i'lo=liblo`i'[1,1]-1.96*sqrt(`=lls`i'[1,1]')
	scalar liblo`i'hi=liblo`i'[1,1]+1.96*sqrt(`=lls`i'[1,1]')
	est restore E`i'
	nlcom _b[Guns:citis]+2*_b[Guns:citisq], post
	matrix libhi`i'=e(b)
	matrix lhs`i'=e(V)
	scalar libhi`i'lo=libhi`i'[1,1]-1.96*sqrt(`=lhs`i'[1,1]')
	scalar libhi`i'hi=libhi`i'[1,1]+1.96*sqrt(`=lhs`i'[1,1]')
}

local i=2000
forvalues i=2000/2002 {
	quietly gsem (Guns -> bminage bmult bssp botheft bdesign bconst btrain, probit) (Guns -> bballist blocks basslt bmags bdtr bmg bss bsbs bsbr baow bgunban bwait bdealer bpriv blic bregis bretent bipcs boci bcci, regress) (Guns -> bipt1, gamma) (blackis urbanis citis citisq south -> Guns, ) if year==`i', latent(Guns) startvalues(iv)
    predict g`i' if e(sample), latent(Guns) se(gs`i')
    est store E`i'
	lincom citis+citisq
    matrix b`i'=e(b)
	matrix c`i'=b`i'[1,59]
	scalar c`i'=c`i'[1,1]
	matrix d`i'=b`i'[1,60]
	scalar d`i'=d`i'[1,1]
	matrix r`i'=b`i'[1,57]
	matrix q`i'=b`i'[1,58]
	scalar r`i'=r`i'[1,1]
	scalar q`i'=q`i'[1,1]
	matrix s`i'=e(V)
	matrix t`i'=s`i'[59,59]
	scalar t`i'=t`i'[1,1]
	scalar c`i'hi=`=c`i''+1.96*sqrt(`=t`i'')
	scalar c`i'lo=`=c`i''-1.96*sqrt(`=t`i'')
	matrix u`i'=s`i'[60,60]
	scalar u`i'=u`i'[1,1]
	scalar d`i'hi=`=d`i''+1.96*sqrt(`=u`i'')
	scalar d`i'lo=`=d`i''-1.96*sqrt(`=u`i'')
	matrix v`i'=s`i'[57,57]
	scalar v`i'=v`i'[1,1]
	matrix w`i'=s`i'[58,58]
	scalar w`i'=w`i'[1,1]
	scalar r`i'hi=`=r`i''+1.96*sqrt(`=v`i'')
	scalar r`i'lo=`=r`i''-1.96*sqrt(`=v`i'')
	scalar q`i'hi=`=q`i''+1.96*sqrt(`=w`i'')
	scalar q`i'lo=`=q`i''-1.96*sqrt(`=w`i'')
	nlcom _b[Guns:citis]-2*_b[Guns:citisq], post
	matrix liblo`i'=e(b)
	matrix lls`i'=e(V)
	scalar liblo`i'lo=liblo`i'[1,1]-1.96*sqrt(`=lls`i'[1,1]')
	scalar liblo`i'hi=liblo`i'[1,1]+1.96*sqrt(`=lls`i'[1,1]')
	est restore E`i'
	nlcom _b[Guns:citis]+2*_b[Guns:citisq], post
	matrix libhi`i'=e(b)
	matrix lhs`i'=e(V)
	scalar libhi`i'lo=libhi`i'[1,1]-1.96*sqrt(`=lhs`i'[1,1]')
	scalar libhi`i'hi=libhi`i'[1,1]+1.96*sqrt(`=lhs`i'[1,1]')
}

local i=2003
forvalues i=2003/2013 {
	quietly gsem (Guns -> bminage bmult bssp botheft bdesign bconst btrain, probit) (Guns -> brifle bballist blocks basslt bmags bdtr bmg bss bsbs bsbr baow bgunban bwait bdealer bpriv blic bregis bretent bipcs boci bcci, regress) (Guns -> bipt1, gamma) (blackis urbanis citis citisq south -> Guns, ) if year==`i', latent(Guns) startvalues(iv)
    predict g`i' if e(sample), latent(Guns) se(gs`i')
    est store E`i'
	lincom citis+citisq
    matrix b`i'=e(b)
	matrix c`i'=b`i'[1,61]
	scalar c`i'=c`i'[1,1]
	matrix d`i'=b`i'[1,62]
	scalar d`i'=d`i'[1,1]
	matrix r`i'=b`i'[1,59]
	matrix q`i'=b`i'[1,60]
	scalar r`i'=r`i'[1,1]
	scalar q`i'=q`i'[1,1]
	matrix s`i'=e(V)
	matrix t`i'=s`i'[61,61]
	scalar t`i'=t`i'[1,1]
	scalar c`i'hi=`=c`i''+1.96*sqrt(`=t`i'')
	scalar c`i'lo=`=c`i''-1.96*sqrt(`=t`i'')
	matrix u`i'=s`i'[62,62]
	scalar u`i'=u`i'[1,1]
	scalar d`i'hi=`=d`i''+1.96*sqrt(`=u`i'')
	scalar d`i'lo=`=d`i''-1.96*sqrt(`=u`i'')
	matrix v`i'=s`i'[59,59]
	scalar v`i'=v`i'[1,1]
	matrix w`i'=s`i'[60,60]
	scalar w`i'=w`i'[1,1]
	scalar r`i'hi=`=r`i''+1.96*sqrt(`=v`i'')
	scalar r`i'lo=`=r`i''-1.96*sqrt(`=v`i'')
	scalar q`i'hi=`=q`i''+1.96*sqrt(`=w`i'')
	scalar q`i'lo=`=q`i''-1.96*sqrt(`=w`i'')
	nlcom _b[Guns:citis]-2*_b[Guns:citisq], post
	matrix liblo`i'=e(b)
	matrix lls`i'=e(V)
	scalar liblo`i'lo=liblo`i'[1,1]-1.96*sqrt(`=lls`i'[1,1]')
	scalar liblo`i'hi=liblo`i'[1,1]+1.96*sqrt(`=lls`i'[1,1]')
	est restore E`i'
	nlcom _b[Guns:citis]+2*_b[Guns:citisq], post
	matrix libhi`i'=e(b)
	matrix lhs`i'=e(V)
	scalar libhi`i'lo=libhi`i'[1,1]-1.96*sqrt(`=lhs`i'[1,1]')
	scalar libhi`i'hi=libhi`i'[1,1]+1.96*sqrt(`=lhs`i'[1,1]')
}

local i=2014
forvalues i=2014/2016 {
	quietly gsem (Guns -> bminage bmult bssp botheft bdesign bconst btrain, probit) (Guns -> brifle bballist blocks basslt bmags bdtr bmg bss bsbs bsbr baow bwait bdealer bpriv blic bregis bretent bipcs boci bcci, regress) (Guns -> bipt1, gamma) (blackis urbanis citis citisq south -> Guns, ) if year==`i', latent(Guns) startvalues(iv)
    predict g`i' if e(sample), latent(Guns) se(gs`i')
    est store E`i'
    matrix b`i'=e(b)
	matrix c`i'=b`i'[1,59]
	scalar c`i'=c`i'[1,1]
	matrix d`i'=b`i'[1,60]
	scalar d`i'=d`i'[1,1]
	matrix r`i'=b`i'[1,57]
	matrix q`i'=b`i'[1,58]
	scalar r`i'=r`i'[1,1]
	scalar q`i'=q`i'[1,1]
	matrix s`i'=e(V)
	matrix t`i'=s`i'[59,59]
	scalar t`i'=t`i'[1,1]
	scalar c`i'hi=`=c`i''+1.96*sqrt(`=t`i'')
	scalar c`i'lo=`=c`i''-1.96*sqrt(`=t`i'')
	matrix u`i'=s`i'[60,60]
	scalar u`i'=u`i'[1,1]
	scalar d`i'hi=`=d`i''+1.96*sqrt(`=u`i'')
	scalar d`i'lo=`=d`i''-1.96*sqrt(`=u`i'')
	matrix v`i'=s`i'[57,57]
	scalar v`i'=v`i'[1,1]
	matrix w`i'=s`i'[58,58]
	scalar w`i'=w`i'[1,1]
	scalar r`i'hi=`=r`i''+1.96*sqrt(`=v`i'')
	scalar r`i'lo=`=r`i''-1.96*sqrt(`=v`i'')
	scalar q`i'hi=`=q`i''+1.96*sqrt(`=w`i'')
	scalar q`i'lo=`=q`i''-1.96*sqrt(`=w`i'')
	nlcom _b[Guns:citis]-2*_b[Guns:citisq], post
	matrix liblo`i'=e(b)
	matrix lls`i'=e(V)
	scalar liblo`i'lo=liblo`i'[1,1]-1.96*sqrt(`=lls`i'[1,1]')
	scalar liblo`i'hi=liblo`i'[1,1]+1.96*sqrt(`=lls`i'[1,1]')
	est restore E`i'
	nlcom _b[Guns:citis]+2*_b[Guns:citisq], post
	matrix libhi`i'=e(b)
	matrix lhs`i'=e(V)
	scalar libhi`i'lo=libhi`i'[1,1]-1.96*sqrt(`=lhs`i'[1,1]')
	scalar libhi`i'hi=libhi`i'[1,1]+1.96*sqrt(`=lhs`i'[1,1]')
}

est restore E1986
test [Guns]blackis = `=r2016'
test ([Guns]citis+2*[Guns]citisq) = `=libhi2016[1,1]'
test [Guns]urbanis = `=q2016'
est restore E2016
test [Guns]blackis = `=r1986'
test ([Guns]citis+2*[Guns]citisq) = `=libhi1986[1,1]'
test [Guns]urbanis = `=q1986'
/*matrix input R=(`=c1986',`=c1986hi',`=c1986lo',`=d1986',`=d1986hi',`=d1986lo' \ `=c1987',`=c1987hi',`=c1987lo',`=d1987',`=d1987hi',`=d1987lo' \ `=c1988',`=c1988hi',`=c1988lo',`=d1988',`=d1988hi',`=d1988lo' \ `=c1989',`=c1989hi',`=c1989lo',`=d1989',`=d1989hi',`=d1989lo' \ `=c1990',`=c1990hi',`=c1990lo',`=d1990',`=d1990hi',`=d1990lo' \ `=c1991',`=c1991hi',`=c1991lo',`=d1991',`=d1991hi',`=d1991lo' \ `=c1992',`=c1992hi',`=c1992lo',`=d1992',`=d1992hi',`=d1992lo' \ `=c1993',`=c1993hi',`=c1993lo',`=d1993',`=d1993hi',`=d1993lo' \ `=c1994',`=c1994hi',`=c1994lo',`=d1994',`=d1994hi',`=d1994lo' \ `=c1995',`=c1995hi',`=c1995lo',`=d1995',`=d1995hi',`=d1995lo' \ `=c1996',`=c1996hi',`=c1996lo',`=d1996',`=d1996hi',`=d1996lo' \ `=c1997',`=c1997hi',`=c1997lo',`=d1997',`=d1997hi',`=d1997lo' \ `=c1998',`=c1998hi',`=c1998lo',`=d1998',`=d1998hi',`=d1998lo' \ `=c1999',`=c1999hi',`=c1999lo',`=d1999',`=d1999hi',`=d1999lo' \ `=c2000',`=c2000hi',`=c2000lo',`=d2000',`=d2000hi',`=d2000lo' \ `=c2001',`=c2001hi',`=c2001lo',`=d2001',`=d2001hi',`=d2001lo' \ `=c2002',`=c2002hi',`=c2002lo',`=d2002',`=d2002hi',`=d2002lo' \ `=c2003',`=c2003hi',`=c2003lo',`=d2003',`=d2003hi',`=d2003lo' \ `=c2004',`=c2004hi',`=c2004lo',`=d2004',`=d2004hi',`=d2004lo' \ `=c2005',`=c2005hi',`=c2005lo',`=d2005',`=d2005hi',`=d2005lo' \ `=c2006',`=c2006hi',`=c2006lo',`=d2006',`=d2006hi',`=d2006lo' \ `=c2007',`=c2007hi',`=c2007lo',`=d2007',`=d2007hi',`=d2007lo' \ `=c2008',`=c2008hi',`=c2008lo',`=d2008',`=d2008hi',`=d2008lo' \ `=c2009',`=c2009hi',`=c2009lo',`=d2009',`=d2009hi',`=d2009lo' \ `=c2010',`=c2010hi',`=c2010lo',`=d2010',`=d2010hi',`=d2010lo' \ `=c2011',`=c2011hi',`=c2011lo',`=d2011',`=d2011hi',`=d2011lo' \ `=c2012',`=c2012hi',`=c2012lo',`=d2012',`=d2012hi',`=d2012lo' \ `=c2013',`=c2013hi',`=c2013lo',`=d2013',`=d2013hi',`=d2013lo' \ `=c2014',`=c2014hi',`=c2014lo',`=d2014',`=d2014hi',`=d2014lo' \ `=c2015',`=c2015hi',`=c2015lo',`=d2015',`=d2015hi',`=d2015lo' \ `=c2016',`=c2016hi',`=c2016lo',`=d2016',`=d2016hi',`=d2016lo')*/
matrix input R=(`=liblo1986[1,1]',`=liblo1986hi',`=liblo1986lo',`=libhi1986[1,1]',`=libhi1986hi',`=libhi1986lo' \ `=liblo1987[1,1]',`=liblo1987hi',`=liblo1987lo',`=libhi1987[1,1]',`=libhi1987hi',`=libhi1987lo' \ `=liblo1988[1,1]',`=liblo1988hi',`=liblo1988lo',`=libhi1988[1,1]',`=libhi1988hi',`=libhi1988lo' \ `=liblo1989[1,1]',`=liblo1989hi',`=liblo1989lo',`=libhi1989[1,1]',`=libhi1989hi',`=libhi1989lo' \ `=liblo1990[1,1]',`=liblo1990hi',`=liblo1990lo',`=libhi1990[1,1]',`=libhi1990hi',`=libhi1990lo' \ `=liblo1991[1,1]',`=liblo1991hi',`=liblo1991lo',`=libhi1991[1,1]',`=libhi1991hi',`=libhi1991lo' \ `=liblo1992[1,1]',`=liblo1992hi',`=liblo1992lo',`=libhi1992[1,1]',`=libhi1992hi',`=libhi1992lo' \ `=liblo1993[1,1]',`=liblo1993hi',`=liblo1993lo',`=libhi1993[1,1]',`=libhi1993hi',`=libhi1993lo' \ `=liblo1994[1,1]',`=liblo1994hi',`=liblo1994lo',`=libhi1994[1,1]',`=libhi1994hi',`=libhi1994lo' \ `=liblo1995[1,1]',`=liblo1995hi',`=liblo1995lo',`=libhi1995[1,1]',`=libhi1995hi',`=libhi1995lo' \ `=liblo1996[1,1]',`=liblo1996hi',`=liblo1996lo',`=libhi1996[1,1]',`=libhi1996hi',`=libhi1996lo' \ `=liblo1997[1,1]',`=liblo1997hi',`=liblo1997lo',`=libhi1997[1,1]',`=libhi1997hi',`=libhi1997lo' \ `=liblo1998[1,1]',`=liblo1998hi',`=liblo1998lo',`=libhi1998[1,1]',`=libhi1998hi',`=libhi1998lo' \ `=liblo1999[1,1]',`=liblo1999hi',`=liblo1999lo',`=libhi1999[1,1]',`=libhi1999hi',`=libhi1999lo' \ `=liblo2000[1,1]',`=liblo2000hi',`=liblo2000lo',`=libhi2000[1,1]',`=libhi2000hi',`=libhi2000lo' \ `=liblo2001[1,1]',`=liblo2001hi',`=liblo2001lo',`=libhi2001[1,1]',`=libhi2001hi',`=libhi2001lo' \ `=liblo2002[1,1]',`=liblo2002hi',`=liblo2002lo',`=libhi2002[1,1]',`=libhi2002hi',`=libhi2002lo' \ `=liblo2003[1,1]',`=liblo2003hi',`=liblo2003lo',`=libhi2003[1,1]',`=libhi2003hi',`=libhi2003lo' \ `=liblo2004[1,1]',`=liblo2004hi',`=liblo2004lo',`=libhi2004[1,1]',`=libhi2004hi',`=libhi2004lo' \ `=liblo2005[1,1]',`=liblo2005hi',`=liblo2005lo',`=libhi2005[1,1]',`=libhi2005hi',`=libhi2005lo' \ `=liblo2006[1,1]',`=liblo2006hi',`=liblo2006lo',`=libhi2006[1,1]',`=libhi2006hi',`=libhi2006lo' \ `=liblo2007[1,1]',`=liblo2007hi',`=liblo2007lo',`=libhi2007[1,1]',`=libhi2007hi',`=libhi2007lo' \ `=liblo2008[1,1]',`=liblo2008hi',`=liblo2008lo',`=libhi2008[1,1]',`=libhi2008hi',`=libhi2008lo' \ `=liblo2009[1,1]',`=liblo2009hi',`=liblo2009lo',`=libhi2009[1,1]',`=libhi2009hi',`=libhi2009lo' \ `=liblo2010[1,1]',`=liblo2010hi',`=liblo2010lo',`=libhi2010[1,1]',`=libhi2010hi',`=libhi2010lo' \ `=liblo2011[1,1]',`=liblo2011hi',`=liblo2011lo',`=libhi2011[1,1]',`=libhi2011hi',`=libhi2011lo' \ `=liblo2012[1,1]',`=liblo2012hi',`=liblo2012lo',`=libhi2012[1,1]',`=libhi2012hi',`=libhi2012lo' \ `=liblo2013[1,1]',`=liblo2013hi',`=liblo2013lo',`=libhi2013[1,1]',`=libhi2013hi',`=libhi2013lo' \ `=liblo2014[1,1]',`=liblo2014hi',`=liblo2014lo',`=libhi2014[1,1]',`=libhi2014hi',`=libhi2014lo' \ `=liblo2015[1,1]',`=liblo2015hi',`=liblo2015lo',`=libhi2015[1,1]',`=libhi2015hi',`=libhi2015lo' \ `=liblo2016[1,1]',`=liblo2016hi',`=liblo2016lo',`=libhi2016[1,1]',`=libhi2016hi',`=libhi2016lo')
matrix input S=(`=r1986',`=r1986hi',`=r1986lo',`=q1986',`=q1986hi',`=q1986lo' \ `=r1987',`=r1987hi',`=r1987lo',`=q1987',`=q1987hi',`=q1987lo' \ `=r1988',`=r1988hi',`=r1988lo',`=q1988',`=q1988hi',`=q1988lo' \ `=r1989',`=r1989hi',`=r1989lo',`=q1989',`=q1989hi',`=q1989lo' \ `=r1990',`=r1990hi',`=r1990lo',`=q1990',`=q1990hi',`=q1990lo' \ `=r1991',`=r1991hi',`=r1991lo',`=q1991',`=q1991hi',`=q1991lo' \ `=r1992',`=r1992hi',`=r1992lo',`=q1992',`=q1992hi',`=q1992lo' \ `=r1993',`=r1993hi',`=r1993lo',`=q1993',`=q1993hi',`=q1993lo' \ `=r1994',`=r1994hi',`=r1994lo',`=q1994',`=q1994hi',`=q1994lo' \ `=r1995',`=r1995hi',`=r1995lo',`=q1995',`=q1995hi',`=q1995lo' \ `=r1996',`=r1996hi',`=r1996lo',`=q1996',`=q1996hi',`=q1996lo' \ `=r1997',`=r1997hi',`=r1997lo',`=q1997',`=q1997hi',`=q1997lo' \ `=r1998',`=r1998hi',`=r1998lo',`=q1998',`=q1998hi',`=q1998lo' \ `=r1999',`=r1999hi',`=r1999lo',`=q1999',`=q1999hi',`=q1999lo' \ `=r2000',`=r2000hi',`=r2000lo',`=q2000',`=q2000hi',`=q2000lo' \ `=r2001',`=r2001hi',`=r2001lo',`=q2001',`=q2001hi',`=q2001lo' \ `=r2002',`=r2002hi',`=r2002lo',`=q2002',`=q2002hi',`=q2002lo' \ `=r2003',`=r2003hi',`=r2003lo',`=q2003',`=q2003hi',`=q2003lo' \ `=r2004',`=r2004hi',`=r2004lo',`=q2004',`=q2004hi',`=q2004lo' \ `=r2005',`=r2005hi',`=r2005lo',`=q2005',`=q2005hi',`=q2005lo' \ `=r2006',`=r2006hi',`=r2006lo',`=q2006',`=q2006hi',`=q2006lo' \ `=r2007',`=r2007hi',`=r2007lo',`=q2007',`=q2007hi',`=q2007lo' \ `=r2008',`=r2008hi',`=r2008lo',`=q2008',`=q2008hi',`=q2008lo' \ `=r2009',`=r2009hi',`=r2009lo',`=q2009',`=q2009hi',`=q2009lo' \ `=r2010',`=r2010hi',`=r2010lo',`=q2010',`=q2010hi',`=q2010lo' \ `=r2011',`=r2011hi',`=r2011lo',`=q2011',`=q2011hi',`=q2011lo' \ `=r2012',`=r2012hi',`=r2012lo',`=q2012',`=q2012hi',`=q2012lo' \ `=r2013',`=r2013hi',`=r2013lo',`=q2013',`=q2013hi',`=q2013lo' \ `=r2014',`=r2014hi',`=r2014lo',`=q2014',`=q2014hi',`=q2014lo' \ `=r2015',`=r2015hi',`=r2015lo',`=q2015',`=q2015hi',`=q2015lo' \ `=r2016',`=r2016hi',`=r2016lo',`=q2016',`=q2016hi',`=q2016lo')
matrix define RS=(R,S)
drop if year<1986
drop if state~="Alabama"
mkmat year, mat(Y)
clear
matrix define out=(RS,Y)
matrix colnames out = citislo citislolo citislohi citishi citishilo citishihi blackis blackishi blackislo urbanis urbanishi urbanislo year 
svmat double out, names(col)
twoway rcap citislohi citislolo year if year>1985 || scatter citislo year if year>1985, ytitle("Marginal Effect of Liberalism at Liberalism=-1")
twoway rcap citishihi citishilo year if year>1985 || scatter citishi year if year>1985, ytitle("Marginal Effect of Liberalism at Liberalism=1")

/*twoway rcap citishi citislo year if year>1985 || scatter citis year if year>1985, ytitle("Liberalism Coefficient")
twoway rcap citisqhi citisqlo year if year>1985 || scatter citisq year if year>1985, ytitle("Liberalism Squared Coefficient")*/
/*figure 6*/
twoway rcap blackishi blackislo year if year>1985 || scatter blackis year if year>1985, ytitle("Percentage Black Coefficient")
twoway rcap urbanishi urbanislo year if year>1985 || scatter urbanis year if year>1985, ytitle("Urban Population Coefficient")

est table E1986 E1990 E1995 E2000 E2005 E2010 E2015 E2016, se stats(bic)
